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Abstract /Summary : 

We  extend  in  this  paper  the  analysis  of  a  posteriori 
estimates  of  the  space  discretization  error  presented  in  a 
previous  paper  [3]  for  time-independent  space  meshes.  In 
the  context  of  the  model  problem  studied  there,  results 
are  given  relating  the  effectivity  of  the  error  estimator 
to  properties  of  the  solution,  space  meshes,  and  manner  in 
which  the  meshes  change.  A  procedure  based  upon  this  theory 
is  presented  for  the  adaptive  construction  of  time-dependent 
meshes.  The  results  of  some  computational  experiments  show 
that  this  procedure  is  practically  very  effective  and  sug¬ 
gest  that  it  can  be  used  to  control  the  space  discretiza¬ 
tion  error  in  more  general  problems. 


1.  Introduction 


The  finite  element  method  of  lines  (FEMOL)  is  the 
numerical  procedure  used  to  approximate  the  solution  of 
time-dependent  partial  differential  equations  (PDEs)-in 
which  a  problem  is  first  reduced,  through  finite  element 
discretization  of  the  space  variables,  to  that  of  solving 
a  system  of  ordinary  differential  equations  (ODEs).  Full 
advantage  can  then  be  taken  of  the  efficient  and  reliable 
ODE  software  available  to  complete  the  computational  pro¬ 
cedure  . 

In  order  to  achieve  a  given  level  of  accuracy  in  the 
FEMOL,  it  is  prudent  to  refine  space  meshes  in  areas  of 
greater  solution  activity.  If  this  activity  is  localized 
in  different  areas  of  space  at  different  points  of  time, 
it  then  is  desirable  to  modify  the  space  mesh  at  various 
times,  and  to  have  decisions  concerning  such  modifications 
made  adaptively  by  the  computer  during  the  solution  process. 
Some  of  the  applications  for  which  this  procedure  is  suit¬ 
able  are  continuum  models  of  enzyme-catalyzed  chemical 
reactions,  the  transmission  of  electrical  nerve  impulses, 
the  propagation  of  flame  fronts  in  combustion,  and  flow  in 
porous  media. 

We  mention  here  that  the  subject  of  mesh  modification 
in  Galerkin  procedures  other  than  the  FEMOL  for  the  solution 
of  time-dependent  PDEs  has  attracted  much  recent  research. 
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The  use  of  finite  elements  both  in  time  and  space  to  modify 
meshes  was  analyzed  in  [11],  and  schemes  implementing  ver¬ 
sions  of  this  technique  can  be  found  in  [5]  and  [6].  The 
concept  of  "moving"  finite  elements  which  deform  continu¬ 
ously  in  time  was  introduced  in  [13],  and  several  interest¬ 
ing  experimental  results  were  presented  in  [10]  and  [14]. 
This  approach  was  recently  analyzed  in  [3],  where  general 
types  of  continuous  and  discontinuous  mesh  modifications 
were  studied,  and  certain  symmetric  a  priori  error  estimates 
were  given. 

This  paper  is  a  continuation  of  a  previous  paper  [3], 
in  which  for  time-independent  space  meshes  we  analyzed  a 
posteriori  estimators  of  the  space  discretization  error 
(i.e.  the  total  error,  provided  that  the  resulting  ODE 
system  is  solved  exactly)  in  the  FEMOL  solution  of  para¬ 
bolic  PDEs .  An  a  posteriori  estimator  based  upon  computable 
local  residuals  was  shown  to  be  effective,  provided  certain 
conditions  on  the  solution  regularity,  mesh  family,  and 
mesh  size  were  satisfied.  Under  less  restrictive  conditions 
a  modification  of  this  principal  estimator  was  shown  to  be 
effective . 

These  results  are  extended  here  to  the  setting  of  the 
FEMOL  in  which  the  space  mesh  is  allowed  to  undergo  discon¬ 
tinuous  transitions  during  the  solution  process.  An  a 
posteriori  estimator  analogous  to  that  given  in  [3,  sec¬ 
tion  5]  is  shown  to  be  effective  in  this  setting,  provided 
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that  regularity  assumptions  on  the  solution  related  to  those 
given  in  [3,  section  5]  hold,  the  discrete  data  at  the  tran¬ 
sition  time  points  are  properly  chosen,  and  the  mesh  transi¬ 
tions  satisfy  certain  conditions. 

The  principal  tool  needed  in  carrying  out  this  extension 
is  the  a  priori  estimate  presented  in  Thm.  3.2.  This  result 
requires  assumptions  related  to,  but  in  some  ways  less  re¬ 
strictive  than  those  given  in  [8].  In  particular,  no  rela¬ 
tions  need  hold  on  the  locations  of  points  in  one  mesh  to 
those  in  the  successive  mesh.  This  freedom  is  desirable 
when  considering  adaptive  schemes  for  the  construction  of 
such  meshes.  The  result  of  Thm.  3.2  is  also  important  in 
that  it  leads  to  apparently  new  optimal  order  a  priori 
error  estimates,  one  of  which  is  demonstrated  in  Thm.  3.3. 

The  a  posteriori  analysis  presented  forms  the  basis 
of  an  adaptive  scheme  given  for  the  control,  as  well  as 
the  estimation  of  the  space  discretization  error.  The 
isolation  of  this  component  of  the  tdtal  error  provides 
insight  as  to  how  such  estimates  and  adaptive  procedures 
can  be  used  with  state  of  the  art  ODE  software  to  estimate 
and  control  the  total  error  in  the  FEMOL. 

All  results  are  given  here  in  the  context  of  the 
model  problem  studied  in  [3]  (cf.  [4]  for  extensions 
mentioned  in  [3]).  The  full  description  of  this  problem, 
necessary  preliminaries,  and  a  discussion  of  some  of  the 
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results  for  changing  meshes  are  given  in  section  2. 
Sections  3  and  4  contain  the  a  priori  and  a  posteriori 
estimates  mentioned  above.  In  section  5  under  conditions 
less  restrictive  than  in  section  4  it  is  shown  that  a 
modification  of  our  principal  estimator,  which  takes  into 
account  the  errors  introduced  by  the  mesh  transitions,  is 
effective.  An  adaptive  scheme  is  presented  in  section  6 
and  applied  to  examples  related  to  those  considered  in 
[3,  section  7], 

In  order  to  facilitate  a  first  reading  of  this  paper, 
some  basic  notions  concerning  a  posteriori  error  estimates 
in  the  setting  of  changing  meshes  are  presented  here,  and 
a  description  of  how  these  estimates  can  be  used  in  an 
adaptive  procedure  to  control  the  space  discretization 
error  is  given.  The  discussion  is  specifically  oriented 
to  the  experiments  described  in  section  6,  which  were 
directed  toward  the  evaluation  of  va: ious  aspects  of  a 
posteriori  estimates  and  adaptive  procedures.  By  reading 
this  description  and  skipping  to  section  6,  one  should 
obtain  a  good  overview  of  the  motivations  behind  the 
somewhat  technical  development  of  the  theory  in  sections 
2-5  and,  more  importantly,  see  how  these  ideas  can  be  used 
in  practice. 

Let  u  =  u(t,x)  be  the  solution  of  the  model  problem 
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(1.1) 


r 


u  =  f-Lu  =  f(t,x)  +  [a ( x ) u  ]  -  b ( x ) u ;  0  <  x  <  1, 

X  X 


u(t,x)  =  g ( t , x )  ; 


u  (  0  ,  x  )  =  u  0  ( x )  ; 


°  <  -t  S  Tfi.ial, 


x  =  0 , 1 ,  0  <  t  5 


0  <  x  <  1  , 


x  FI MAI 


where  a  >  0,  b  5  0,  u^,  g,  and  f  are  given  functions. 

For  simplicity  we  assume  here  that  these  functions  are  both 
smooth  and  compatible,  but  the  theory  does  not  require  this. 

The  FEMOL  approximation  U  of  u  is  the  solution  of 
a  weak,  or  integral  form  of  eqs.  (1.1),  and  is  defined  through 
the  partitioning  of  the  time  and  space  intervals 


( 0  ,  T 


FINAL 


)  and  (0,1).  The  procedure  used  to  create  these 


partitions  and  obtain  U  is  based  upon  the  attempt  tc 

control  for  each  t  i  [0,TrTMAI  1 

FINAL 


(1.2)  III  e  ( t ,  •  )  HI  5  (  [ 

W  0 


a(x)e^(t,x)dx^'/“ 


where  e  =  u  -  U,  according  to  the  following  principle. 
For  any  given  TOL  >  0,  the  partitions  of  (0,1)  used  in 
the  FEMOL  should  be  chosen  or  modified  at  various  times  so 
that 


(1.3)  III  e  ( t ,  •  )|||  ~  TOL;  Vt  €  CO  »TFINAL]  . 

Required  input  for  the  procedure  are  partitions 


0 


(1.4)  A, 


fn  0  0 
{0  =  XQ  <  < 


<  X, 


-  1} 


and 


(1.5)  A 


{0  =  T0  <  T1  < 


<  T  =  T  > 

M+l  1  FINAL 1 


The  points  {T  }  n  are  times  when  discontinuous  mesh 

m  m=  0 , M+l 

changes  are  permitted,  although  they  need  not  occur.  The 
first  step  in  the  algorithm  is  the  determination  of  a 
continuous,  piecewise  linear  function  U(C,*)  on  A^ , 
which  is  chosen  to  approximate  u(0,*)  =  u ^ ( • ) .  The  norm 
III  e  (  0  ,  •  )||f  is  then  estimated  by 


(1.5) 


r  No  9^1/2 

E<o)  ■  U  ”«i}  ’ 


where  the  computable  local  error  indicators  {r)n  .}._ 


are  defined  by 


0j Jj=l,N 


0 


(1.7)  T) 


Oj 


|x?-x? 

Lj _ izL 


fX"  n+x? 

12al  L-l-1 


1 


0 

xi-l 


(Lu-LU) 4 (0 ,x)dx;  1  £  j  £  N, 


Based  upon  the  size  of  each  a  decision  of  whether 


to  accept 


s0 


'03 

and  U(0,*)  or  to  modify 


Aq  and  recompute 


U(0,*)  is  made  by  the  computer.  Having  obtained 


and 


U(0,O,  the  FEMOL  solution  U  =  U(t,x)  is  then  uniquely 
determined  as  a  function  which  is  continuous  and  piecewise 


/ 


linear  in 

x  and 

smooth  in 

t  on  each  subinterval 

in  the 

following 

manner . 

For 

-  €  [Ta 

,TX],  V(t 

,  • )  is  characterized  by  the 

vector 

of  its  nodal  values  {U ( t , x° )  }  .  _  v  - ,  which  is  the 

3  1  “-1-  1  '»Q 

solutior.  of  an  Nq  -  1  dimensional  ODE  initial  value  problem. 
In  order  to  isolate  the  space  discretization  error  and  means 
of  estimating  it,  very  small  error  tolerances  were  used  in 
a  time  discretization  scheme  in  the  computational  examples 
so  that  these  ODEs  were  essentially  solved  exactly.  The 
norm  jjj  e  (T^  ,  •  )l||  is  then  estimated  by 


(1.3) 


E<V 


V/: 


to; 


where  N ^  =  N  ^ =  { C  =  xjj  <  x^  <  .  .  .  <  x?j  =  1}  =  Aq, 

and  the  indicators  {p.  .  }.  ,  ..  are  defined  by 

-‘■3  3  ~  1 2. 


(1.9)  r^. 


1  1  i 2  xt 


x • -x .  . 

3  3-1 


1  x  1 

|X .  , +x . 

12a/-^- 1 


3 


1 

X3-l 


(bV  +  LU-f )"(T1  ,x)dx; 

t- 


1  £  3  £  N. 


Usine  {n,  a  decision  of  whether  or  not  to 

6  '!]  ]=1  ,N 


accept  A, 


and  U(T^,*)  is  then  made.  If  A^  is  modified, 


new  initial  or  transition  data  U(T^,*)  and  U^CT^,*),  which 
are  continuous  and  piecewise  linear  on  the  new  mesh  A^ , 
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are  determined.  This  process  is  effected  through  the  solution 

of  linear  algebraic  equations.  The  error  is  again  estimated 

with  E(T^),  defined  now  on  the  new  mesh  a  ,  and  the 

modification  process  is  permitted  to  occur  once  more. 

For  each  m  =  1,M  the  process  of  solving  for 

t  €  [T  ,T  ,,  ]  an  '!  -  1  dimensional  ODE  oroblem  with  initial 
m  m+I  m 

data  {'J(T  jX1?)}.  ,  ,  is  likewise  carried  out,  and  the 

m  j  j=l,W  -1 
J  J  m 

decisions  of  whether  or  not  to  modify  the  partition 


.  r .  _  m  m  m 

( 1 .  x  0 )  A  =  { 0  =  xn  <  x,  <  .  .  .  <  x„ 

m  0  1  h 


m 


are  similarly  made  by  the  computer  after  checking  the  error 

a  ■P  n  i 

> 

'm 


indicators  {q  .  } .  ,  '  ,  defined  by 

mj  1=1,  N 


m 


(1-n>  ^m  j 


I  m  m  1 2  x . 
X  .  -X  .  ,  r  1 

3  1-1  '  J 


12a 


m  ,  m 
X  .  ,  +x  . 

_ 1 


m 


(U  +LU-f)2(T  ,x)dx ; 
_  m 


■j-1 


1  S  j  <  N  , 
J  m 


and  the  estimator  E(Tm)  °f  III e C ,  *  )!!j  ,  given  by 


(1.12) 


:(T  )  =  fNf  q2.V/2. 

m  Vj=l  mV 


The  quantity  HI  e  ( t ,  ■  )l|l  is  clearly  defined  for  all 

t  €  C 0  .Tt,,...  T  ]  and,  in  an  obvious  way,  the  estimator  E(‘) 
FINaL 

can  be  extended  to  a  piecewise  smooth  function  on 


& 


g 

In  order  to  assess  the  performance  of  E(-)  and  the  adapt!’ 
procedure^  we  define  the  effectivity  ratio  G(-)  bv 

(1.13)  G  ( t )  =  E  ( t )  /  ET-  ( t )  ;  7 1  €  [0,TrTV  ,  T],  where 

(1.14)  eTRtj ( t )  5  ;'e(t,  •)'!!. 

The  process  of  estimating  :|| e (  •  ,  *  )l||  with  E(-)  on 

[0,TpTj^T]  is  effective  if  l/0(t)  is  uniformly  bounded 

on  [0,Tr-.,A.],  and  it  is  desirable  that  3(t)  would 

approach  1  for  each  t,  were  the  procedure  repeated  with 

successively  smaller  input  tolerances  TOL.  This  in  turn 

would  lead  to  an  effective  adaptive  procedure,  since  the 

decisions  at  each  T  are  made  in  an  attempt  to  control 

m 

E(T  )  so  that  E(t)  ~  TCL  for  each  t  €  [C,TrTMA,]. 
m  FINAL 

The  theoretical  results  in  sections  u  and  5  related  to 
the  effectivity  of  the  estimator  E(-)  are  given  for 
generally  constructed  meshes  and  are  asymptotic  in  nature. 

In  section  6  experiments  testing  the  applicability  of  these 
conclusions  in  the  above  adaptive  procedure  are  presented 


and  discussed. 
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2 .  Preliminaries,  Problem  Description,  and  Discussion 
Throughout  this  paper  the  convention  that  the  variables 
i,  j,  t  ,  m,  M,  n,  N,  p,  and  a  take  only  integer 
values  shall  be  used.  We  begin  by  introducing  certain  nota¬ 
tion  and  concepts  needed  for  the  study  of  the  model  problem. 
For  more  details  we  refer  to  [3]  and  [12]. 

Let  I  =  (0,1)  c  o  5  T4  <  T  <  » ,  the  open  inter¬ 

val  J  c  (0,T),  real  s  >  0,  and  V  be  an  arbitrary  Hilbert 
space . 

We  denote  by  C^( J ; V) {H^ ( J;V) }  the  space  of  functions 
which  are  continuous  {square  integrable}  on  J  with  range 
in  V,  and  by  HS(J;V)  the  corresponding  Sobolev  space  of 
order  s.  When  J  =  I  and  V  =  R2'  we  suppress  J  and  V 

g 

and  denote  by  H  the  usual  Sobolev  space  of  order  s  on  I 

with  norm  ||  •  ||  ,  and  by  the  subset  of  functions  in  H1 

s  u 

which  vanish  at  0  and  1.  H2s,S(J)  =  H°(J;H2s)  fl  HS(J;H°) 
is  the  nonisotropic  Sobolev  space  of  order  s  on  J  x  I  with 
norm  II  '  H^s  ,s  * 

Throughout  this  paper  a  and  (3  denote  real  numbers 
related  to  the  data  of  the  model  problem  and  satisfy 

(2.1)  -1  <  p  <  a;  3  <  a  S4  ;  p,a  t  integer  +  1/2. 

For  sufficiently  smooth  functions  a  and  b  on  I 
satisfying 

(2.2)  0<asa(x)<a<«} 


(2.3) 


0  <  b(x)  <  b  < 


Vx  e  i. 
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we  consider  the  model  problem:  Find  u  =  u(  ;u^,f)  satisfy¬ 
ing 


(2.4) 


~  (t,x)  +  Lu(t,x)  =  f(t,x);  t,x  €  (0,T)  *  I, 

O  L 

i 

\u(t,0)  =  C  =  u(t,l);  t  €  ( 0 , T) , 

^u(0,x)  =  uQ (x) ;  x  €  I, 


where  L  is  the  positive  self-adjoint  operator  defined  by 


(2.5)  Lu  =  - ( au  )  +  bu 

X  X 

with  domain  D(L)  =  (v  £  H2  :  v(x) I  n  ,  =  0). 

x=0,l 


More  generally,  the  domains  of  fractional  powers  of  L  are 
Hilbert  spaces  defined  by 

(2.6)  P(L(P+1)/2)  =  (v  €  H?+1  :LPv(x)Ix=q  1  =  0 

Vp  satisfying  0  s  p  <  p/2  +  1/4}. 

The  data  and  f  in  eas.  (2.4)  are  required  to 

satisfy 

(2.7)  U(),f  €  Da,P(0,T)  =  { v ,  g  e  HP  +  1  x  Ha,a/2(0,T)  : 

M^R(v,g)  =  0}, 


q 


where  with  £  ' 

=  o 

for  q  <  p 

and  [u] 

denoting 

the 

3  =P 

,.CR 

is  defined 

integral  part  of 

any 

real  number 

a, 

by 

fD 

(2.8)  Mp  ' v,g)  = 

I 

[0/2+1/4] 

Xv(x)+  l 

’(-l)PLi-p 

3P-- 

at^(0’x)l 

r  il 

x=0 , 1 

i  =  0 

P  =  1 
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Setting 

(2.9)  Ha,i3(T,8,T)  =  {u  €  Ha  +  2’(a  +  2)/2(Tft,T)  :  u  =  u(  ;v,g)  solves 

eas.  (2.4)  for  some  data  v,g  €  Ea,^(0,T)}, 

CL  8 

we  define  for  u(  ;v,g)  €  H  ,P(T.,{,T) 


( 0  , T) 


(2.10)  Ca(P,T*,T,u)  = 


|v|,p  +  l  +  l|g|la,a/2  ’ 


(1+T-(a+l)/2v„CR 


if  T*=0, 

if  T*  >  0  . 


We  have  the  following  useful  result. 


Theorem  2 . 1  (cf.  [3,  Thm.  3.4]).  Let  uQ,f  €  Da’^(0,T), 
with  p  =  a  if  =  0,  and  set  T*  =  0  if  p  =  a.  Then  eas. 

(2.4)  have  a  unique  solution  u(  ;uQ,f)  €  Ha’^(T.,.,T)  and  3 
a  positive  constant  C  such  that 


u 


(T*,T) 

a+2 , (a+2)/2 


S  CC  (P,T.,T,u) 


Using  the  restriction  (2.1)  in  [3,  Thms .  2. 1-2. 3],  we 
also  have 


Theorem  2.2.  Let  u  €  Ha’^(Tft,T)  and  JC(T.,{,T). 
Then  3  a  positive  constant  C  such  that  for 


(2.11)  j  =  0,  1,  and  2, 


a^u 

atJ 


€  H 


and 


d2u 


"it? 


(a  +  2-  2j ) , (a  +  2-2j  )/2 


(«2-2j),(«2-2j)/2fJ)nH0(J.Hl) 

-  C||ulla*2,(C,t2)/2’  and 


(2.12)  for  j  =  0  and  1,  €  C  (J;H' 


J  n 


and  sup 
t€J 


E  ~§T~J  u  ( t ,  •  ) 

J  ot  ‘a  +  l-: 


ct  +  2  ,  (a  +  2)/2  • 


As  in  [3],  in  order  to  demonstrate  the  principal  a 
posteriori  error  estimates  in  section  4  we  require  the  use  of 

q  Q 

subclasses  of  H  ’  (T,.,T),  consisting  of  functions  satisfying 
one  or  both  of  the  following  additional  properties. 

a  8 

(2.13)  For  5^  >  0  we  say  that  u  €  H  ,H(T.,5,T)  satisfies 

property  1(6^)  if  for  each  t  6  [Tft,TJ,  llu^(t ,  •  ) !!  >  0 


sup  \ 
t  €  [T...  >T]  '• 


C  (P  , T..j , T, u) 


Fu  (t , 

1  XX 


.  5  V 


(2.14)  For  6 ^  >  0  we  say  that  u  €  H  ’ H(T...  ,T)  satisfies 
property  11(6  )  if  for  each  t  €  [T..,,T]  3  an  open 
interval  I  c  I  such  that 

gT(t,u)  s  inf  { lu  (t,x) ! •  |l.  |1/2}  >  0  and 

L  X€T  xx  t 


gTJ(t.u) 

t€[TAP,T]  gL(t,u)  5  62 


,  where 


gu(t,u') 


sup(|u  (t,x)[  +  lu  ( t  ,  X )  f  } 
^  XX  ’  XXX 


For  5  >  0,  6  >  0  we  then  set 
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(2.15)  tf“,P(T.,,T)  =  {u  €  tfa,P(T,.,T)  :u  satisfies  I(6n)},  and 
°1  "  "  1 


(2.16)  ( T ,  T )  =  {u  €  H*’P(TA,T)  :u  satisfies  IK  60  )  > 

0-|  ^  o  2  **  l. 


Our  numerical  analysis  is  based  upon  the  following  weak 
form  of  eqs.  (2.4).  Find  u  :  (0,T)  -*•  satisfying 

(2.17)  <u  (t)  ,*>  +  B(u(t),<*>)  =  <f(t)  ,<#>>*,  Vt,<$  €  (0,T)  x  Hq  , 


(2.18) 


u(0)  =  uQ  , 


where  <•,•>  denotes  the  inner  product  on  I  and  the  sym¬ 

metric  bilinear  form  B  is  defined  by 


(2.19)  B(u,v)  =  <au  ,v  >  +  <bu,v>; 

X  X 


Vu , v  e  H 


By  (2.2)  and  (2.3)  3  positive  constants  C^,  such  that 


(2.20)  C  2 1|  w  ||  ^ 


l|w||^  =  B( w , w)  <  C1!|w||^;  Vw  €  Hj! 


Let  T  and  P  denote  the  families  of  mesh  partitions  of 
(0,T)  and  I,  respectively.  For 


{0 


T0  <  T1  <  •■•  <  Tm+1  =  T}  5  T\  M  =  M(A  )  >  1,  and 


(a  } 


m 


m 


=0 ,M( AT) 


,  with 


V  =  <°  =  xo  <  X1  < 


•  •  <  xJJ  =  1}  €  P;  m  f  [0,M],  Nm  >  2, 


m 
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we  write 


Jm 

=  <VW> 

n 

B 

It) 

CTm’W- 

CJ 

m 

=  CVW’ 

Jm  s 

"m'Vl1' 

T 

m 

=  T  .  -  T  , 
m+1  m 

> 

,  m  m  s 

- 

!Tj  i 

m  m 

xj  "  Xj-1’ 

h 

— m 

>  i  ..m  j 

=  mm  |  I  •  |  , 

1<  j  <N  3 

m 

hm 

=  max  |  Z? | 

l<j<Nm  J 

m  / 

0  <  m  <  M(  a  )  > 


1  5  3  :  Nm, 


and 


h  =  h(A)  =  max  T  h  . 

0<m<M(A  ) 

T  T 

For  A  and  A  as  above,  we  denote  by  a  *  A  the  parti- 

M(AT) 

tion  of  (Q,T)  x  I  into  the  7  N  two  dimensional  subin- 

_n  rn 

m-  0 

tervals  of  the  form 


{J  x  O  T  ’ 

3  m=0,M(AT) 

j=l,N 
J  m 

In  general,  A  %  A  for  p  \  q ,  but  by  setting  A  =  A„  for 
&  pq  J  °  m  0 

T 

m  =  1,M(a  ),  our  results  will  directly  reduce  to  those  shown 

in  [3]  for  time-independent  space  meshes. 

We  say  that  the  family  T  *  P  of  such  Cartesian  product 

partitions  of  (0, T)  *  I  is  y-regular ,  if  3y  i  0  and 

T 

>  0  such  that  for  each  A  *  A  €  T  *  P 
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(2.21)  T  2  A,h2r;  0  <  m  <  H(AT), 

m  1  m 

and  that  T  x  P  is  <,y-regular  if  T  *  P  is  y-regular  and 

3<  >  1,  A0  >  0,  and  A^  >  0  such  that  for  each 
T 

A  *  A  €  T  x  p 

(2.22)  h?  <  A2hq;  0  <  p  <  q  <  M(a?), 

and 

(2.23)  h  >  A,hK;  0  s  m  <  M(aT). 

— m  j  m 


(2.24)  Unless  otherwise  specified,  throughout  the  remainder 
of  this  paper  C  shall  denote  a  positive  generic 
constant  which  depends  in  general  on  the  functions  a 


and  b,  the  constants 

V 

A  2  > 

A3>  and  T,  but  is 

independent  of  M, 

m, 

Tm  ’ 

A, 

m 

h  ,  T* ,  and  the 

functions  Uq,  f, 

and 

u . 

For  each  A  €  P,  S(A 
m 

)  c 
m 

is  the 

subset  of  functions 

restrictions  to  any 

I? 

3 

are 

linear . 

The  operators  P™, 

P‘£,  and  L&  mapping  H^,  Hg ,  and  S(Am),  respectively, 

to  S(a  )  are  defined  by 
m 


5(  v  ,  ) 


17 


La  is  a  positive  self-adjoint  ODerator  with  domain  S(A  ) 
m  *  m 

and  satisfies 


(2.26)  L.  (P%)  =  p"J(Lv); 

m  1  u 


Vv  €  P(L) , 


(2.27)  j|  L“ i«!lQ  S  C II ^  If  ^  S  fU  I!  Q 
m  m 


V*  €  S ( A  )  . 

m 


We  have  the  following  approximation  results. 


Lemma  2 . 1  (of.  [3,  Lemmas  2.1,  2.3]).  There  exists  a  posi¬ 
tive  constant  C  such  that  for  each  A  6  P 

m 


(2.28) 


II  (I-P  )vil 
1  s 


cCS!|v!l  ; 
m  cr 


Vv  €  V(La/2), 
0SSS1SOS2, 


(2.29) 


(I-P™)vi1  5  Ch®U„  , 

0  0  m  cr 


Vv  e  v ( lc/ 2 ) , 


0  5  a  S  2,  o  t  1/2, 


(2.30)  5  ch’llvll^,  hMjb) 


L  =  j+o/2,  j  2  0,  1  S  a  S  2  . 


For  each  A  *  A  €  T  x  P ,  we  define  the  FEMOL  solution 


M<  AT) 

U  :  [ 0 ,T  )  -  U  S ( A  ) 
m=0  m 

by 

(2.31)  U(t)  =  (At);  t  €  cJm  for  some  m  €  [0,M], 


where 
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U°  :  JQ  -  S< A0> 


is  defined  by 


C2.32)  <U°(t),*>  +  B(U°(t),*) 


<f(t),$>;  Vt  ,<j>  €  Jg  xStijj), 


(2.33) 


and 


U°(0) 


u. 


[piV 

if 

r  -a 

o  o 
c 
o 

if 

U311  :  J_  -►  S (a  )  for  1  s  m  <  M(a^) 
m  m 

is  defined  by 

(2.34)  <U™(t)  ,<*>  +  B(Um(t)  ,*)  =  <f(t),*>;  U,*  (J^xS(4m), 


(2.35) 


Um(T)  =  L“1P^(f(T)  -U*_1(Tm)) 

m  Am  0  m  t  m 


From  the  definition  of  L  ,  we  see  that  eq.  (2.35)  is 


equivalent  to 


m 


B(Um(T  ),*)  =  <f(T  )  -U™-1(T  )  ,<t>>\  V*  €  S  ( A  ), 

m  m  t  m  m 


and  it  therefore  follows  from  eqs .  (2.32),  (2.34)  and  (2.35) 
that 


(2.36) 


1  im  t/?(t) 

+ 

t-T 

m 


pmUm_lfT  )  • 
0ut  ' 


1  <  m  <  M(  A  ) 


Eqs.  ( 2 . 31) - ( 2 . 3 5 )  are  equivalent  to  a  sequence  of 
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uniquely  solvable,  iteratively  defined  ODE  initial  value 

problems.  At  each  transition  time  T  ,  the  transition  data 

m 

(new  initial  data)  is  determined  from  the  solution  of  an  N  -  1 

m 

dimensional  matrix  equation. 

Remarks .  In  the  case  of  pure  mesh  derefinement ,  i.e. 


A  c  A  ,  for  some  m  >  1,  then  S(a  )  c  S(A  -,  )  and  we  have 
m  m-1  m  m-1 


from  eqs 

.  (2.34) 

and  (2.35)  that 

(2.37) 

Um( T  )  = 

m 

P™Um~1(T  ). 

1  m 

In 

the  case 

of  pure  mesh 

refinement , 

i.e.  A 

m 

m-1 

then  S ( a  )  o  S(a  ,  )  and  it 
m  m-1 

follows  that 

(2.38) 

pin  lyin^  }  . 

1  m 

and 

lim  U®(t)  = 

t-T+  t 

U™-1(T  ). 
t  m 

im 

For 

u  €  Ha 

,P(T*,T)  and 

T 

each  A  x  A 

€  T  x  p 

we  define 

(2.39) 

P1(t) 

:  Pm- 
*1* 

t  6  °J 

m 

for  some 

m  €  [0 ,M] , 

(2.40) 

e(t) 

=  u ( t )  -  U(t) 

9 

t  €  [0 ,T)  , 

(2.41) 

p(t) 

-  P^(t)u(t)  - 

u(t)  ; 

t  €  [  T  ,  T )  , 

(2.42) 

9  ( t ) 

=  P  (t)u(t)  - 

U(t)  ; 

t  €  [T*,T). 

We 

shall  use  the  following  notational 

convention.  For 

1  <  m  <  M, 
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In  section  3  it  will  be  shown  with  y  <  1  and  certain 
a  =  a(y),  that  if  u  €  Ha,^(TA,T)  and  T  x  ?  is  a  <,y- 
regular  mesh  family,  then  the  choice  (2.35)  of  discrete  transi¬ 
tion  data  yields  an  a  priori  estimate  of  the  form 

(2.44)  ||  0(t )  ||  +  ||e  (t)||  =  0(h°)  as  h  -  0 ; 

L  U  ITl  HI 


Vt  €  CT.,{  ,T  5  eacn  m  €  [0, K],  and  some  cr(y,a)  >  1.  The 
result  (2.44)  and  the  restriction  (2.23)  for  certain  k  will 
lead  to  statements  in  section  4  about  the  effectivity  of 
estimating  IMIjr  on  CTft,T)  with  an  a  posteriori  estimator 

E(  •  )  . 

We  note  that  in  the  FEMOL,  the  y-regularity  restriction 
(2.21)  on  T  x  p  is  apparently  necessary,  not  only  for  effec¬ 
tive  error  estimation  but  also  convergence.  We  refer  to  the 

counterexample  in  [8,  section  4]  in  which  a  sequence  of  meshes 

4 

satisfying  t  =  o(h  )  as  h  -*■  0  led  to  convergence  to  the 
mm  m 

wrong  solution  of  the  heat  equation  in  E''- . 

We  next  remark  that  in  order  for  the  estimator  E(-)  to 
perform  well  on  [T.,{,T),  T  x  p  must  satisfy  a  restriction 
like  (2.22).  If,  for  example,  a  sequence  of  partitions  of 
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(0,T)  x  I  was  used  for  which  A*  =  {0  <  T/2  <  T}  and  A  ^ 
remained  fixed,  while  h^  -+  0,  then  for  T/2  5  t  <  T  we 
would  have  E(t)  =  O(h^)  as  h,  -»  0,  while  'i  e  ( t ) I]  >  C  >  0 
in  general . 

The  restrictions  (2.22)  and  (2.23)  would  not  be  needed  if 
our  only  interest  vjere  in  obtaining  a  priori  error  estimates. 
If  T  *  P  is  y-regular  with  y  5  1,  u  €  Ha,^(T... ,!)  for 
certain  a  =  a(y),  and  the  discrete  transition  data  choice 
(2.35)  is  made,  then  in  section  3  we  shall  show  estimates  of 
the  form 

f  !|  e(  t )  ||  „  =  0(h)  "l 

(2.45)  <  2-v  [  as  h  -+  0 ;  Vt  6  [T,.  ,T). 

Ule(t)H0  =  0(h2  Y)J 

More  generally,  by  replacing  {S(Am)}m=Q  ^  in  the  FEMOL  with 
{ ( A  )}  n  where  for  i  >  1  S1(A  )  c  H^. 

denotes  a  closed  subspace  possessing  the  approximation  pro¬ 
perties  of  the  piecewise  i-degree  polynomials  on  A^,  it 

can  be  shown  (cf.  [4])  that 

r  II  e  ( t )  ||  E  =  0  (h1 ) 

<2.46)  \  i+l-r 

Ule(t)||0  =  0  (h1  1  r) 

We  note  that  an  estimate  of  the  form 

( ft  2  11/2  i 

(2.47)  ||  e  ( t )  ||  0  +  M  ||  e  ( s )  ||  ^  dsj  =  0(hx)  as  h  -  0 


r 

► 


follows  from  a  recent  result  given  in  [3],  where  with  a  dif¬ 
ferent  choice  of  transition  data  and  a  1-regular  family 
T  x  p  it  was  also  assumed  that  {S1(i  ,  )  1  S~(A  )}  n  M 

has  certain  approximation  properties,  which  is  a  restriction 
on  the  locations  of  points  in  any  mesh  to  those  of  points  in 
the  successive  mesh. 

For  computational  reasons,  in  the  adaptive  scheme  pro- 

oosed  in  section  6  it  is  assumed  that  each  partition  A  is 

m 

obtained  from  the  previous  partition  A^  .  through  a  combina¬ 
tion  of  refinement  and  derefinement,  but  the  theory  permits 
for  instance  meshes  which  satisfy 


(2.48) 


m 

<  x^  < 


m-1 
X1  < 


<  x 


m 

M 


m 


vm_1  T 

•1  XN  , -1 
m-1 


in  which  case  S(A  ,)  0  S  C  A  )  =  {0}. 

m-1  m 


3.  A  Priori  Estimates 


Let  u(  ;un,f)  he  the  solution  of  eqs.  (2.4),  for  each 

T 

A^  x  a  8  T  x  P  let  U  be  the  FEMOL  solution  defined  by  ecs. 
( 2 . 31 ) - ( 2 . 3 5  )  ,  and  e,p,  and  8  be  as  defined  in  (2.40)  — 
(2.42).  The  following  theorem  is  a  direct  consequence  of  the 
result  in  [3 ,Thm. 4 . 1 ] . 

Theorem  3.1.  Let  u(  ;Ug,f)  8  Ha  ’  “  (T,.  ,T)  ,  where  a 

and  3  satisfy  (2.1).  Then  3  a  positive  constant  C  such 
T 

that  for  each  A  *  A  8  T  x  P  with  T^  >  T.,.  and  all 
t  8  CJ0  n  [T.,s  ,T) 

(3.1)  j|  6  ( t )  |l  E  +  l|  e  t  ( t )  j|  q  2  ChJ“2Ca(  ?,TiV,T,u). 

By  using  Thm.  3.1,  the  transition  data  choice  (2.35), 

and  the  embedding  results  in  Thm.  2.2  on  each  subinterval 

J  c  (Tju.T),  we  will  exploit  the  local  time  character  of 
m  •• 

the  error  and  be  able  to  show 

Theorem  3.2.  Let  u(  ;uQ,f)  8  Ha  ’  p  ( T... ,  T )  ,  where  a 
and  8  satisfy  (2.1),  and  T  x  p  be  a  <,y-reeuiar 
family.  Then  3  a  positive  constant  C  such  that  for  each 
AT  x  a  8  T  x  P  with  T^  >  T... ,  each  m  8  [0,M],  and  ail 
t  8  cJm  n  [T.v  ,T ) 

(3.2)  ||  0<t )  |!e  +  II  e  t  ( t )  I!  0  <  Ch^“2~YCa(P,TA,T,u). 


Proof.  Since  y  ?  0, 


the  result  (3.2)  for  m 
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follows  directly  from  Thm.  3.1,  and  it  therefore  suffices  to 
demonstrate  (3.2)  for  arbitrary  n  satisfying  1  £  m  *5  M. 


Suppose  we  have  shown  that  positive  constants 


and 


c  - 


m 


a-2, 


(3.3)  II  e*.(  t  )l!  g  -  l|et(Tm)!lC  +  Clnm  !l  u!1  a  + 2  ,  (a  +  2  )  /  2  ’  and 


a-2. 


J 


(3.4)  ||9  (t)ll  _  <  Vl!et(Ti),|0*Clh^l'ulhh,(«42)/2>- 


We  then  claim  that  for 

(3.5)  j[  e  t  ( t )  ||  Q  <  |let(Tp|l0 

(3.6)  !|0  (t)||  £  <  CQ{  ||  et  (T“ )!) 


all  t  €  J 


m 


m  „  7  J  a 

o-2. 


+  C1  u'*a  +  2  ,  (a  +  2  )/2 


,  and 


m  „  J. 
.  +  C,  l  h,  Hull 

0  1  •  n  J 

u  3  =  1 


a  +  2, (a  +  2  )/ 2 


}. 


If  m  =  1  there  is  nothing  to  show,  as  (3.5)  and  (3.6) 
are  then  just  (3.3)  and  (3.4).  If  m  >  2  we  can  estimate 
the  first  terms  on  the  right  hand  sides  of  (3.3)  and  (3.4) 
by  following  with  m  replaced  by  m  -  1  the  same  steps 
which  led  to  (3.3),  since  by  the  nctational  convention  (2.43) 

(3.7)  II  e  (T~)  II  n  =  lim  II  e .  ( t )  ||  n . 

t  m  u  t  u 

1m 


Repeating  this  argument  a  total  of  m  -  1  times  yields  (3.5) 
and  (3.6). 


As  eq.  (3.7)  also  holds  for  m  =  1,  we  have  by  Thm .  3.1 
that  3  a  positive  constant  C  such  that 


(3.8)  l|et(Tl)!l0  5  Ch0~2_Y  ca<3’T*’T’u)- 

Combining  (3.5),  (3.6),  and  (3.8),  it  follov/s  that  3  a 
positive  constant  C  (independent  of  m)  such  that  for 
all  t  €  CJ 

m 


(3.9) 


Meet)!!. 


r  o  m 

:(t)!l0  £  C^hJ-2-YCa(?,Tft,T,u)+ 


a-2-y 


h?Y/r.) 
3  I/ 


1/2  1  /  9  J- 

t.  2||u||  2 
3  l|u|la  +  2,(a+2)/2 


Using  m  £  M  and  the  « , y-regular itv  conditions  (2.21) 
and  (2.22),  (3.9)  becomes 


(3.10)  |l  ©  ( t )  ||  _  +  II  e  ( t )  ||  n  £  Ch  °"2"r(c  (  P ,  T ... ,  T ,  u ) 

ti  T  u  m 


M  ,  J, 
+  I  x!  Dull 


3  =  1 


a+  2 , ( a+  2  )  / 2 


} 


by  the  M-dimensional  Cauchy-Schwar z  inequality 
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<  Ch 


a-2-y 


m 


M 

U  )  +  (  l  T  .  ) 

j=l  ^ 


1/2 


(  .^  (  ■u'!a+2,(a+2)/2)  ) 


Ch 


a-  2  -  y 


m 


f  1/2 

<Ca(^,Tft,T,u)  +  (T-T1 ^  ■  i|  u  I'  ‘ 


a+2, (a+2)/2 


} 


£  Chr2"MCa(e’T-T'u)  *  ^u'h  +  2  ,  ( a+2  )  /  2 


{■ 


(T,.,T) 


by  Thm.  2.1 


<  Ch^_2~YCa(p,T)V,T,u), 


which  is  the  desired  result  (3.2). 

We  shall  now  demonstrate  (3.3)  and  (3.4),  first  for 

t  =  T  and  then  for  t  €  J  .  We  have  bv  (2.43)  and 
m  m 

observation  (2.36)  that 


(3.11)  e.(T  ) 
t  m 


u  (T  ) 
t  m 


P0Ut'1(Tm) 


=  P™  lim  (u  ( t)-U  ( t) )  +  (I-P^)u.(T  ) 
u  t  t  u  t  m 

m 

=  P^e  (T-)  +  (I-P^)u(T  ). 

0  t  m  0  t  m  . 

Hence,  by  the  definition  of  P  ™ :  H°  -*■  S(A  )  and  the  triangle 
inequality 


(3.12)  lIVV'lo  s  l'et(Tm),!0  +  »(I-pS)lit(Tm)!l 


by  (2.1)  and  (2.29) 
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s  llet^T~)  IIq  *  Ch“-2  SU£  llut(t)lla.2 

t  J 

m 


by  (2.12) 


s  +  Chm'2 


J 


m 


a+  2  ,  (  a+  2  )  /  2  ‘ 


3y  (2.43),  (2.17),  (2.35),  and  (3.11)  we  have  that 


(3.13)  0 (T  ) 

m 


•L”1P'V(T  ) 

Am  0  t  m 


-L ^pjje  (T“) 
A  0  t  m 
m 


Therefore  by  (2.27)  and  the  definition  of  P^ 

(3.14)  ||  G  (T  )  ||  _  S  C||  e  (T~ )  II  n  . 

m  h  t  m  0 

Now,  it  follows  from  the  t-independence  of  the  operator  L 
and  subtraction  of  eq.  (2.34)  from  eq.  (2.17)  that  for 
j  =  0  or  1 

/aj+l  \  /aj+i  \ 

(3.15)  s),*>  *  ! 

»s,*  €  Jm  «  S(4m). 

With  <t>  =  20s(s)  and  j  =  0,  eq.  (3.15)  becomes 

(3.16)  2 1|  0g  ( s )  ||  q  +  ~  ||  0  ( s )  ||  ^  =  2<^ps(s)  ,0s(s)^>  . 

For  arbitrary  t  €  J  ,  integration  of  eq.  (3.16)  on  (Tm,t) 
and  application  of  the  Cauchy-Schwarz  inequality  yield 
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(3.17)  ||9(t)||2  s  I]  9  (T  )||  l  +  2 

£  m  E 


I"  »<PT-I)us(s)1'o-  "es's)»0ds 


m 


r  .i« 


(  s )  l|  2ds 
s  0 


m 


!i  ( P^-I  )us  (  S  )  |)  gds  . 


m 


Hence,  3  a  positive  constant  C  such  that 

m  J 

(3.18)  l|0(t)||E  <  II  0  (Tm)  ||  E  +  C||  (P^-I)utH0^0 

by  (2.1),  (2.11),  and  (2.30)  <  ||0(Tm)|lE  +  Ch°~2 

Now,  by  the  definitions  (2.39)  -  (2.42) 


!iuii 


j 

m 

a+2 , ( a+  2 ) / 2  ‘ 


(3.19)  e . (t )  =  -P.(t)  +  9  .  ( t )  ;  Vt  6  J  . 

t  t  t  m 

Setting  <t>  =  29g(s)  in  eq.  (3.15)  with  j  =  1,  we  get 

(3.20)  |rlles<s)ll20  *  2i|9s(s)ll|  =  2<pss(s),6s(s>). 

As  in  (3.17)  and  (3.18),  integration  of  eq.  (3.20)  on 

(T  ,t)  for  t  €  J  yields 
m  m 


(3.21)  II et ( t ) II 0  <  l|et(Tm)||0  +  C||(p^-I)uttll0”c 


S  «VV"o  *  Chr2MaT2,Ca*2)/2 
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By  (3.19),  (3.21),  and  the  triangle  inequality  we  have  for 
t  €  J 

m 

(3.22)  tl  e  t  ( t )  l|  0  <  ||et(Tm)||Q  +  ||pt(t)||0  +  ilot(Tm>  ||  Q 

+  Chm  a+2  ,  (  a+2  )  /  2 


^  e 


t'V", 


+  2  sun  ||  ( I-P?)u  ( s )  || 

S  €  J  1  S  U 

m 


+  Ch 


a- 2 
‘m 


Hull 


J 

m 

a+2 , (a+2 )/2 


by  (2.1),  (2.12),  and  (2.28) 

<  l|et(Tm),|0  *  Ch“-2||un^2j(a(J)/2. 

Combining  (3.12),  (3.14),  (3.18),  and  (3.22)  completes  the 
proof  of  (3.3),  (3.4),  and  hence,  Thm.  3.2  . 

We  now  show  the  a  priori  estimates  mentioned  in  section  2 . 

Theorem  3.3.  Let  u(  ;Ug,f)  €  Ha ’ ^ (TA ,T) ,  where  a 
and  8  satisfy  (2.1),  and  T  *  P  be  a  y-regular  family 
with  y  S  a  -  3.  Then  3  a  positive  constant  C  such  that 
for  each  A ^  *  A  €  T  x  P  with  T^  >  T.,{  and  all  t  €  [TjV,T) 

(3.23)  f| e ( t ) || E  <  Ch(A)CaO,Trt,T,u) 


and 
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(3.24)  ||  e  C  t )  ||  0  <  Cha-2~Y(A)Ca(3,TA,T,u). 

Proof .  Let  t  €  CJm  H  [T...  ,T)  for  any  m  €  [0,M].  Then 

(3.25)  ||  e  ( t )  ||  E  =  li-p(t)  +  0(t)||E  <  ||  (P^-I)u(t)  ||E  +  I)  6  ( t )  ||  £ 

by  (2.12)  and  (2.28)  <  Ch(A)  sup  l|u(  s )  |l  +  ||  9  ( t )  !| 

s  €  C  T  ,.c ,  T  )  2  E 

by  (2.1),  (2.12),  and  Thm.  2.1  <  Ch(A)C  (3,T,.,T,u)  +  ||0(t)l|-. 

a  t 

Similarly , 

(3.26)  ||  e  ( t )  ||  g  <  Ch2(A)Ca(p,T,,,T,u)  +  !!  ©  ( t )  (|  Q 

by  (2.27)  <  Ch2(A)Ca(  3,T*,T,a)  +  C||e(t)!!r 

Examination  of  the  proof  of  Thm.  3.2  shows  that  the 

k , Y-regularity  condition  (2.23)  was  not  used,  and  that 

(2.22)  was  only  used  in  obtaining  (3.10)  from  (3.9).  Therefore, 

replacing  h  with  h(A)  in  (3.10)  yields 
m 

(3.27)  ||  0  ( t )  ||  r  <  Cha_2_Y(A)C  (  P,TA,T,u). 

lj  a 

Since  y  S  a-  3  by  assumption,  combining  (3.25),  (3.26), 
and  (3.27)  yields  the  desired  results  (3.23)  and  (3.24). 

Remarks .  If  a  =  4,  the  results  (3.23)  and  (3.24) 


become 
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(3.28)  h"(1~Y)||e(t)||0  +  ||e(t)|iE  £  ChC^C  e,T*  ,T,u) ;  Vt  €  CTa,T), 

and  if  p  =  4  and  T}.{  =  0  (smooth  and  compatible  data), 

(3.23)  and  (3.24)  become 

(3.29)  h_(1"Y)||e(t)||0  +  ||e(t)  |1E  £  Ch||u||^’T) 

£  Ch{||u0H5  +  !|f||^’T)};  Vt  6  [0  ,T )  . 

The  estimate  of  Hell-  is  of  optimal  order  in  h,  while  the 
estimate  of  ||e||g  is  only  °f  optimal  order  when  the  number 
of  mesh  changes  is  finite  (y=0),  a  result  given  in  [7]  with 
a  different  choice  of  data  at  the  transition  times. 
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4 .  A  Posteriori  Error  Estimates 
We  begin  this  section  by  extending  the  notions  of  upper 
and  asymptot ically  exact  estimators  presented  in  [3]  to  the 
setting  of  time-dependent  space  meshes. 

For  A^  *  A  €  T  *  P  with  ,  let 


(4.1) 

.  T 

g:  A 

X 

A  -*• 

A 

g 

a 

f  m . 
=  {g  } 

T 

9 

and 

m=0 ,M( A 

) 

(4.2) 

T 

<t>  :  A 

X 

A  -*■ 

4 >A 

A 

,  m , 

=  U  } 

T 

5 

with 

A 

m=  0 ,M( A 

) 

(4.3) 

U 

g  m  € 

c° 

(CJ 

m 

n 

[T.V,T)  ; 

H° )  ^ 

A 

c° 

1 

j 

> 

m  €  [ 0 ,M] , 

(4.4) 

4>  m  € 

(CJ 

m 

n 

CT*,T); 

S(A  ))  J 

m  J 

$  €  C°([T*,T);  Hq) . 

We  say  that  a  quantity  E(*,g,<£)  is  an  upper  estimator 
for  on  [I,..  ,T)  if  3  positive  constants  ,  h* 

such  that 

(4.5)  for  each  A^  *  A  €  T  x  P  with  >  T.v  and  h(A)  5  C^h* , 

(4.6)  ||#(t)-0  A  ( t )  ||E  <  C2  E(t,gA,0A);  Vt  ?  [Ta  ,T)  , 

where  and  are  as  in  (2.24)  and  h*  depends  in  general 

on  the  solution  class  in  which  $  resides  (i.e.g.  on 
a,  p,  6^,62)  and  on  the  family  T  x  P  (i.e.g.  on  k,y)- 

If  under  the  same  hypotheses  3  a  number  a  >  0  such 


that 
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(4.7)  ||$(t)-<pA(t)||E  =  E(t,gA,4>A)(l  +  0(ha))  as  h  -  0; 

Vt  €  [ T ,  T )  , 

where  the  constant  in  0(ha)  is  as  in  (2.24),  then  we  say 
that  E(*,g  ,<f>)  is  an  asymptotically  exact  estimator  for 
ll$-*llE  on  [T*,T). 

Let  u(  ;Ug,f)  be  the  solution  of  eqs.  (2.4)  and  for 
T 

each  A  x  A  €  T  x  P  let  U  denote  the  FEHOL  solution 

defined  by  eqs.  ( 2 . 31 ) - ( 2 . 3 5 ) .  We  now  define  our  principal 

a  posteriori  estimator  for  j|u-U||  . 

For  each  m  €  [0,M]  and  t  €  cJm  Cl  [T*,T),  the  local 

indicators  {n  are  defined  by 

-  m3  1=1 ,Nm 


(4.8)  r1^j(t,f-Ut,U)  = 


I^l2 


12a 


X™ 


.  ,  +X3? 
3-1  3 


m(Ut+LU-f)  ( t ,x)dx , 

Ij 


and  the  associated  estimator  E  by 


f  m  n 

(4.9)  E(t ,f-U. ,U)  =  {  l  rV\(t,f-U 

r  kj=l 

E(*,f-U^,U)  is  a  piecewise  smooth  function  on  [T.V,T) 
is  based  upon  local  computable  residuals,  and 


(4.10)  E(t,f-Ut,U) 


0(h(  A )  )  as  h(  A)  ->  0  ;  Vt  €  [T*,T) 


34 


Under  suitable  hypotheses,  by  applying  Thm.  3.2  and 

results  shown  in  [3]  for  time-independent  space  meshes  to 

c 

each  of  the  subintervals  J  f!  [T...  ,T),  we  shall  be  able 

m 

to  show  that  E(*,f-U  ,U)  is  an  upper  or  asymptotically 
exact  estimator  for  '] u— U || „  on  [T...  ,T).  The  hypotheses 
leading  to  these  results  are  among  the  following: 


(4.11) 

u 

€  Ha,?(T*,T), 

where  a  and 

p  satisfy 

(2.1), 

(4.12) 

u 

€  ^  c  ’  ^  ( T ,  T )  ; 

51 

5^  >  0  and 

T  x  P 

is  K 

, Y-regular , 

with  T,  >  T,. 
1 

for 

each 

AT  X  A  6  T  X 

(4.13) 

u 

6  ( T ... ,  T )  ; 

51’52  " 

61  >  0,  52 

>  o, 

(4.14) 

< 

+  Y<4-2  (=> 

Y  €  [0,1)  and 

K  € 

[1,2)  by  (2.1)). 

As  in 

[3],  the  results  are  demonstrated 

by  comparing 

l|e||E 

to 

!l  P II E  and  E(-, 

f-U  ,U)  to 

E( • , f- 

ut,P1u) 

,  where 

E( ‘ ,f- 

ut,P 

^u)  is  defined 

.  as  in  (4.9), 

with 

Ut(t) 

and 

U(t) 

replaced  by  u^(t) 

and  ?^ ( t )u( t ) . 

We  have  the  following  sequence  of  lemmas. 


Lemma  4 . 1  (cf.  [3, Lemma  5.1]).  There  exists  a  positive 

T 

constant  C  such  that  for  each  A  x  a  €  T  x  p  and 

t  €  CJ  0  [ T a , T ) ;  any  m  €  C 0 , M] 
m  ’• 


|E(t,f-Ut,U)  -  E(t,f-ut,P1u)  |  <  Chm  { |1  ©  ( t  >  It  £  +  !|et(t)||0>. 
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Lemma  4 . 2  (cf.  C 3 , Lemma  5.2]).  Assume  C u  .11)  .  Then 
E  (  •  ,  f-u^  ,  P^u)  is  an  upper  estimator  for  jj  o  |!  on  [T.V,T), 
with  h*  =  1. 

Lemma  4 . 3  (cf.  [3, Lemma  5.3]).  Assume  (4.11),  (4.12), 

and  the  result  of  Lemma  4.2.  Then  3  a  positive  constant 

C  such  that  for  each  A^  x  A  €  T  x  P  and  t  €  CJ  D  [T...  ,T) 

m 

any  m  €  [0,M] 

ChK 

E(t,f-ut,P1u)  ?  C  a  (  P  ,  T ... ,  T ,  u )  . 

Lemma  4 . 4  (cf.  [ 3 , Lemma  5.4]).  Assume  ( 4 . 11 ) - ( 4 . 14 ) . 
Then  E  ( • , f-u^ , Pn u )  is  an  asymptotically  exact  estimator 
for  i|p||  on  [T,..,T),  with  a  =  1  -  k/2  and  h*  = 

The  following  results  can  now  be  shown. 

Theorem  4.1.  Assume  (4.11),  (4.12),  and  (4.14).  Then 
E(*,f-U^,U)  is  an  upper  estimator  for  Hu-U'lj,  on  [T.,.,T), 
with  h*  =  6  ^  ^ ,  where  |i  =  (a-2-K-y)/2. 

Theorem  4.2.  Assume  ( 4 . 11 ) - (4 . 14 ) .  Then  E(*,f-U  ,U) 
an  asymptotically  exact  estimator  for  !|u-U!|r  on  [T...,T), 
with  a  -  |i  and  n-  =  minlS^  ’52  where 


is 


p.  =  ( 0-2-K-y)  / 2  . 


36 


Proof  of  Thms.  4.1  and  4.2.  We  follow  closely  the  steps 


used  in 

[3  , Thms . 

5.1,  5.2]. 

Assume  (4 

.11 )  ,  (4.12),  and 

(4.14), 

T 

let  A 

x  A  €  T  x  P , 

and  let 

t  6  cj  n  [  T ... ,  T ) 

m 

for 

some  m 

€  [  0  ,  M  ]  . 

To  simplify  notation 

we  suppress  t 

and 

also  write 


(4.15) 


E]_  =  E(t,f-Ut,U), 

<  E2  =  E(t,f-ut,P1u), 
^Ca  =  Ca(p,T;,,T,u). 


By  Lemma  4.3  we  may  write 


(4.16)  1 


Ei-E2 


¥ 


E1  1E1'E2 
<  1  <  1  +  ■  lrl~ 

t2  2 


9 


which  after  application  of  Lemma  4.1  and,  again,  Lemma  4.3 
becomes 


(4.17)  1 


C6,  h1-< 
1  m 


{||e|lE+! 


to 


}  5 


*1 


<  1  + 


C61  h1_< 
1  m 


.{||ellE+i|eti 


By  (4.11),  (4.12),  and  (4.14)  we  may  apply  the  result  of 
Thm.  3.2  to  (4.17).  In  view  of  the  definition  of  n,  we  get 


(4.18)  1 


C6.  h2ti+1 
1  m 


1  +  C51h 


2  p.+l 
m 
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Now,  from  the  definition  of  o,  we  may  apply  the  triangle 
inequality  to  e  =  -p  +  0  and  get 


(4.19) 


Up! 


s  I  ei 


s  1 1  Pi 


which  by  Lemma  4.3  and  Thm.  3.2  becomes 


||pliE  HpIL  2 

(,.20)  q  -zf  s  lleilj.  s  ej^  *  C^2m 

Setting  h*  =  we  see  in  (4.18)  that  3  a  positive 

constant  C,  such  that  for  h  2  h  5  C,  h* 

1  ml 

(4.21)  E2  =  E  ( 1  +  0(h^+1 ) )  as  h  -  0, 

and  therefore  from  (4. 20  that 


(4.22)  ||  e  ||  E 


E  (1  +  O(h^) ) 

i  c2 


as  h  -*■  0  . 


By  Lemma  4.2,  the  result  of  Thm.  4.1  follows  immediately 
from  (4.22).  By  assuming  (4.13)  we  have  Lemma  4.4,  which 
applied  to  (4.22)  yields  the  result  of  Thm.  4.2,  since 
0<  n  <  1  -  k/2  <1/2. 
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5.  A  Posteriori  Error  Estimates  Revisited 
In  this  section  we  show  how  the  estimator  E(*,f-U  ,U) 
can  be  modified  to  take  more  fully  into  account  the  evolution¬ 
ary  nature  of  eqs .  (2.4)  and  the  errors  introduced  by  the 
discontinuities  in  the  data  at  the  transtion  time  points. 

Under  less  restrictive  conditions  than  (4.H)  -  (4.14),  this 
modification  is  shown  to  be  an  upper  estimator  for  j|u-U|L 
on  [0 ,T)  . 

We  assume  that 


(5.1)  uQ  €  D(L)  and 

(5.2)  f,  ft  €  C°( [0 ,T) ;  H°), 

from  which  it  follows  that  u(;Ug,f)  €  C^([0,T);  D(L))  and 
u  €  C° ( [ 0 ,T ) ;  H° ) .  Define 

2 

(5.3)  c  =  arr2(a  as  in  (2.2))  <  inf  ^ . 

"  _  <P*0  Hcpllj 

T  c 

For  each  A  *  A  €  T*  P  and  t  €  any  m  €  [0,M], 

we  define 

4 

(5.4)  E(t,f-U,rJ)  =  l  E .  ( t ,  A  ><A),  where 

r  i  =  l  1 

E1(t ,  AT  x  a)  =  E(t,f  -  Ut,U) , 


(5.5) 
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By  assumptions  (5.1),  (5.2),  and  the  definitions  (2 . 31)-( 2. 35) , 
and  (2.43)  we  see  that  E(t,f-U  ,U)  is  well-defined  and 
computable,  given  U. 

We  show  the  following  result. 


Theorem  5.1.  Assume  (5.1),  (5.2),  and  that  T  *  P  is 
y-regular  with  r  S  1.  Then  E(*,f-Ut,U)  is  an  upper 
estimator  for  ||u  -  Ull^.  on  C0,T),  with  h'{  = 


1. 
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Proof.  We  must  show  that  there  exists  a  positive  con- 

T 

stant  C  such  that  for  each  A  x  a  €  T  *  P  and  all 

t  $  CJ  ;  any  m  €  [0,M] 
m 

(5.10)  ||u(t)  -  U(t)|L  5  C  E(t,f  -  U  ,U) . 

t  t 

As  in  [3,  Thm.  6.1],  to  this  end  we  decompose  the  error 
e  =  u  -  U  as 

(5.11)  e  =  e^  +  e^,  where 

(5.12)  e1  =  v  -  U, 

(5.13)  e^=  u-v, 


and  for  each  t  €  [0,T),  v(t)  €  V(L)  is  defined  as  the 
unique  solution  of 


(5.14) 


Lv( t )  =  (f  -  Ut) (t) . 


By  (5.1),  (5.2),  (2.43),  and  the  t-independence  of  the 
operator  L,  we  have  that  v  satisfies 

(5.15)  Lv  (t)  =  (f  -U  Hr);  Vt  6  [0.T),  and 

*■*  U  L 

(5.16)  |-£jV(t)  e  ( cJm; P( L) ) ;  for  j  =  0  and  1  and  all 

m  €  [ 0 ,M ] . 


From  eqs .  (2.17),  (2,18),  ( 2 . 31 ) - ( 2 . 3 5 )  ,  and  the  definitions 
(2.39)  and  (2.43)  we  have  that 
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(5.17)  P1(t)  =  l^-jUCt);  Vt  €  [  0  ,T )  ,  j  =  0  and  1. 

It  follows  directly  from  results  shown  in  [3,  Thm.  6.1]  that 
3  a  positive  constant  C  such  that 

(5.18)  ||e2(T“)!!r, 

and  more  generally  that 

(5.19)  the  desired  result  (5.10)  holds  when  m  =  0. 

In  view  of  the  definitions  (5.4)  —  (5.8)  ,  (5.11),  and  the 

observations  (5.18)  and  (5.19),  the  demonstration  of  (5.10) 

T 

will  be  completed  by  showing  that  for  each  A  x  A  €  T  x  P 

and  t  €  CJ  ;  1  5  m  <  M 

m 

(5.20)  ||e1(t)||E  2  C  E(t,f-Ut,U),  and 

,  ~C(  t-T,  )  /2  9  rn 

(5.21)  lie  (t)ilE  <  C{e  1  ||  e^  (T~ )  il£  +  E2(t,A  x  A) 

+  E^tjA^x  A)}. 

In  view  of  the  result  in  Cl,  Thm.  7.2],  as  applied  in 
[3,  (6.22)],  we  have  the  estimate  (5.20)  by  (5.2),  (5.12), 
(5.14),  and  (5.17)  with  j  =  0. 

In  order  to  show  (5.21), we  subtract  eqs.  (  2.34),  (  2.35) 
from  eqs.  (2.17),  (2.18)  and  get  from  (5.12),  (5.13),  and 
(2.43)  that  for  any  m  satisfying  1  s  m  5  M, 


lim  ||e(t)ILsC{E.(T,,AixA)  +E,  (T,  ,  A 1  x  A  )  }  , 
T  E  2  1  4  1 

1 

t€J„ 


(5.22)  <e2(s),o>  +  B(e2(s),.?)  = 


-<e*(s),o>;  Vs,p  €  Jm  x  H^, 


(5.23) 


e2(T  )  =  e2(T")  +  v(T~)  -  v(T  ) . 

ra  m  m  m 


It  follows  from  (5.9),  (5.14),  (2.43),  and  (2.36)  that 


(5.24)  B(v(T~)  -  v(T  ),<p) 

m  m 


Vo  e  hq, 

0  :  Vo  €  S(A  ), 

’  m 


and  so  by  the  approximation  result  (2.28)  and  a  duality  argu¬ 
ment  there  exists  a  positive  constant  C  such  that 

(5.25)  llv<T->  -  v(Tm)ll0  i  Chml|v(T;)-v(Tm>!lE 

S  cVlIJJI  . 

m  m  U 

From  (5.23),  (5.25),  and  the  triangle  inequality  it  therefore 
follows  that 

(5.26)  He2  (TOIL  <  !ie2(T")|L  +  ChJ|JJn,  and 

m  E  m  E  m  m  0 

(5.27)  il  e 2  ( T  )  li  n  2  ile2(T")l!n  +  Ch2f|J  ||  . 

mu  m  0  m  m  0 


We  now  note  from  (5.17)  with  j  =  1  that 


B(ei(s),o)  =  0 ;  Vs,o  6  J  x  S ( A  ) 
s  m  m 


Hence  ,  it  follows  by  using  a  duality  argument  and  the  defini¬ 
tion  of  E(*,f  -U  ,U  ),  as  the  definition  of  E(*,f-U  ,U) 
s  ss  s  s 

was  used  in  showing  (5.20),  that  for  some  positive  constant  C 
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(5.28) 


«>>»0  s  ChmE(s'fs-Uss’Us>;  Vs  '  V 


As  in  [3,  (6.23),  (6.24)],  we  get  from  (5.28)  and  energy 
arguments  in  eqs .  (5.22),  (5.23)  that 

(5.29)  l|e2(t)!|r  S 


-c( t-T)/ 2  1/?  2 

;  :  *min(2  lie  (T  )||_, 

m  h 


(t-T  ) 
m 


<W 


+  Ch  e 
m 


-c( t-T_ )/2  (ft  -c(Tm-s)^2 


m 


E  ( s , f  -U  ,U  )ds 
’  s  ss  s 


'm 


■) 


1/2 


Vt  €  J  . 
m 


It  can  similarly  be  shown  that 


(5.30)  ||  e  v  t)  II  Q  < 


-c(  t-T  )  /  2  /  ,  (ft  -c(T  -s) 


m 


E  <s'fs-Us5 


,us)dSj1/2); 


Vt  €  CJ 


m 


by  (5.27) 

-$(t-T  )/2  /  „  2  f -S(Tm"s) 

s  e  "  ■('le2<Tm)!l0*Ch2i|Jmll0tChm<  e 


m 


E  (s'fs-Uss 


,us)ds}1/2); 


Vt  €  CJ 


m. 


Combining  (5.26)  and  (5.29),  we  get 
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(5.31)  !|e2(t)HE  5 


-c( t-T  )/2 
m 


21/2lU2<T;)llE*chini|Jm!l0*Chm. 

E2(s,f  -U  ,U  )ds\  1/2  \ 
s  ss  s  /  /’ 


/•t  -c(T  -s) 
m  ■ 


m 


Vt  €  J  . 
m 


In  view  of  the  definitions  (5.6)  and  (5.7),  this  is  the  de¬ 
sired  result  (5.21)  if  m  =  1.  If  m  >  2,  we  can  estimate 
||e2(Tm)||E  in  (5.31)  by  applying  (5.27)  and  (5.29)  with  m 
and  t  replaced  by  m  -  1  and  Tm-  Using  the  y -regularity 
of  the  family  T  x  p  in  this  argument,  we  get 


(5.32)  ||  ez  ( t )  ||. 


,  . 0  -c( t-T  , ) / 2  0 

(2/A  )1/2  •  e  m_1  •  ||eZ(T"  )  || 

1  m- 1  u 


-c( t-T  ) / 2 


+  Ce 


m 


0  -c ( T  -T  ,)/2 

( h  II J  l|  +h2  ,e  ra  ra-1  || J  JL) 
m  m  0  m-i  m-1  0 


-c(t-T  )/2  m 


+  Ce 


m 


•  l  h 

j=m-l  -1 


-c(T  -s) 
m 


J  .  0  (  0  ,  t ) 
3 


E2 ( s , f  -U  ,U  )ds> 17  2 ;  Vt  €  CJ  , 
s  ss  s  /  m 


where  is  as  in  the  y-regularity  definition  (2.21). 

Noting  from  the  definition  (5.3)  of  c  that 


!!  e2  (T~  )  II  Q  <  (l/S1/2)1!e2(T1)'lE, 


(5.33) 
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we  see  that  the  desired  result  (5.21)  follows  from  (5.32)  if 

m  =  2.  If  m  5  3,  we  can  estimate  !!e2(T_  .  )  !|  n  in  (5.32) 

m-1  0 

by  applying  (5.30)  with  n  and  t  replaced  by  m  -  2  and 
T^_ ^ .  Repeating  this  argument  a  total  of  m  -  2  times,  we 
get  from  (5.32)  and  (5.33)  that 


(5.34) 


!!e  ( t )  i|  E  s 


'  1/2  c(t-T.)/2  * 

(2/cA1)i/^-e  1  •  lie2  (T~)  il£ 


+  Ce 


-c( t-T  ) / 2 
m 


m-1  _  -c(T  -T  .  ) / 2 

(h  II J  j|_  +  l  h2e  m  ]  !|J.  iL) 

m  m  0  -  ■>  ->  n 


3=1 


3  'O’ 


-c( t-T  ) / 2  m 
+  Ce  m  •  I  h. 

j=l  : 


■c ( T  -s) 
m 


J.n(o,t) 


E2( s , f  -U  ,U  )ds) 1/2;Vt  €  CJ  , 
s  ss  s  j  m 

The  result  (5.21)  for  m  2  3  follows  immediately  from  (5.34), 

thus  completing  the  proof  of  Thm.  5.1. 

We  now  give  conditions,  which  along  with  (5.1)  and  (5.2) 

are  sufficient  to  show  that  up  to  higher  order  terms  in  h, 

E ( • , f  -  U  , U )  ~  E ( • , f  -  U  ,U)  . 
t  t 


(5.35) 


T,.  €  ( 0  ,T )  and  C*  >  0. 


(5.36) 


I  u  ||  +  l|f(0)||n  +  sup  If  (t)ll  5  C“ . 

0  2  0  t€[0,T]  *  0 


(5.37)  T  x  P  is  a  r-regular  family  with  y  €  (0,1). 


(5.38)  For  each  AixA€TxP,  T  >  and  h(A)  ST4 


1/  2  y 
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The  data  and  f  and  family  T  *  P  are  such  that  for 

each  x  A  €  T  x  P 

f  T 

(5.39)  sup  iU  ( t )  !  +  (j  :|U  (t)!j^dt)1/2 

t€  [T...  ,T)  t  h  J0 

f  i1  9  1/9  it 

+  (  |  !i  U  ( t )  j!  fdt ) 1  SC. 

j  Q  tt  0 

We  have 

Theorem  5.2.  Assume  (5.1),  (5.2),  and  ( 5 . 3 5 ) - ( 5 . 39 ) .  Then 

(5.40)  E(t,f-U  ,U)  =  E(t,f-U  ,U)  +  0(h2_r)  as  h  -*■  0 ;  Vt  €  [T...  ,T)  , 

where  the  constant  in  the  0-term  depends  on  C",  T,  in 

the  y-regularity  definition  (2.21),  and  the  functions  a 
and  b . 

The  proof  of  Thm.  5.2  shall  be  omitted,  as  it  easily 
follows  by  showing 

(5.41)  E(0, Ut(0)  +  LuQ  -  f(0)  ,0)  <  Ch  { :i  uQ  !|  ^  +  '|  f  (  0 )  ||  Q  }  , 

(5.42)  !iJnll0  <  Ch!|Ut(T")!l£;  1  <  m  <  M, 

(5.43)  E(s,f  -U  ,U  )  <  Ch{  j  f  (s)-U  (s)*|  +I|U  (s)|l_>-.  Vs  €  [0,T), 

S  SS  S  S  SS  US  tj 

and  applying  the  Cauchy-Schwarz  inequality  and  y-regularity 

T 

of  T  x  P  to  bound  the  sums  appearing  in  E^CtjA  *  A)  and 

( t ,  A 1  x  a)  . 

T  T 

Suppose  that  for  A  x  a  €  T  x  p,  M  =  M(A  )  2  3  and 

t  €  CJ„.  Since  the  terms  in  E_,  E,,  and  E  have  time 

M  2  3  4 


decay  rates  which  are  pessimistically  small 


any  given 


for 

problem,  in  view  of  Thm.  5.2  one  practical  implementation  of 

E(t,f-U  ,U)  would  be  to  use  only  E(t,f-U  ,U),  and  to 
t  t 

monitor  the  quantities 


(5.44)  E.( t,A  *  A) 


M 


l  h. 
j=M-p  ^ 


-c(t-s) 


j .  n  ( o ,  t ) 
3 


E 2 ( s , f  -U  ,U  )ds} 1/2 
s  ss  s  ' 


and 


(5.45)  E3(t,A  x  1) 


-c( t-T  ) / 2 

e  m  hm!|J  I|n+  r  \fe 

m  m  0  ^ 


M-l  .  -c ( t-T  .  )  / 2 
2  3 


j  =  M-  p  3 


»V«’ 


where  p  =  0,  1,  or  2. 


Also,  since  E9  depends 

of  the  changes  in  time  of  the 

{n  . ( t , f  -  U .  ,U)  } ,  defined  in 
'm  j  t 

be  effectively  monitored  with 
putational  work. 


essentially  on  the  past  history 
local  error  indicators 
(4.8),  both  E^  and  E^  can 
a  minimal  amount  of  extra  com- 
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6 .  Adaptive  Procedure  and  Computational  Examples 

In  this  section  we  outline  many  aspects  of  a  general 
FEMOL  procedure  in  which  decisions  based  upon  a  posteriori 
error  estimates  are  used  to  adaptively  construct  and  modify 
the  space  meshes.  This  procedure  is  contained  in  the  re¬ 
search  program  FEMOL 1,  which  was  written  and  implemented 
on  the  IBM  370  System  in  the  Division  of  Computer  Research 
and  Technology,  NIH.  For  a  discussion  of  the  general  capa¬ 
bilities  and  specific  computational  aspects  of  FEM0L1,  we 
refer  to  [3]  and  [4].  Most  comments  on  the  procedures  em¬ 
ployed  in  FEMOL 1  shall  here  be  restricted  to  those  related 
to  the  adaptivity  in  our  experiments.  These  experiments 
were  designed  for  the  evaluation  of  conclusions  concerning 
a  posteriori  error  estimation  which  are  based  upon  the 
asymptotic  analysis  presented  in  sections  3-5  for  general 
changing  meshes.  The  discussion  is  primarily  oriented  to 
selected  examples,  in  which  the  exact  solutions  are  those  in 
the  examples  presented  in  [3,  section  7],  The  forms  of  thes 
solutions  are  typical  of  those  arising  often  with  systems  of 
linear  and  nonlinear  parabolic  equations  used  in  various 
fields  of  application. 

Before  beginning,  we  mention  that  in  [9]  an  adaptive 
procedure  in  some  ways  similar  to  our  own  was  implemented  fo 
the  FEMOL  solution  of  parabolic  PDEs  in  one  and  two  space 
dimensions.  Using  small,  uniform  time  steps  in  the  backward 


Euler  method  to  integrate  the  resulting  ODE  systems,  the 
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flexibility  of  the  procedure  was  compared  to  that  using  uniform, 
time- independent  meshes  in  a  few  examples. 

The  model  problem  considered  here  is  that  given  in  eqs. 
(1.1),  where  the  functions  a,  b,  and  the  exact  solution  u 
will  be  specified  in  each  of  the  examples,  and  u^,  f,  and 
g  are  chosen  to  be  the  corresponding  smooth,  compatible 
data.  While  the  theory  in  sections  2-5  was  presented  for  the 
case  when  the  boundary  data  g  s  0,  it  holds  more  generally 
for  eqs .  (1.1). 

The  notation  used  here  is  that  of  section  1.  For 

further  clarity  we  recall  that  the  norm  .'!e(t,*)!!l  of  the 

space  discretization  error  e  =  u-U  was  defined  in  (1.2), 

the  initial  mesh  An,  error  indicators  {i..}.,  ..  and 

estimator  ECO)  in  (1.4),  (1.6),  (1.7),  and  each  subsequent 

mesh  A  ,  corresponding  error  indicators  {q  .}.  ,  and 

m’  ^  &  m]  j  =  1 ,N 

J  J  m 

estimator  E(Tm)  used  in  the  adaptive  procedure  in  (1.10), 

(1.11),  (1.12). 


The  evaluation  of  E(*)  and  the  adaptive  procedure  is 
based  upon  examining  the  effectivity  ratio  0( • )  and  the 
true  error  •  )  ,  defined  in  (1.13)  and  (1.14).  While 

the  theory  in  sections  4-5  was  presented  for  the  estimation 
of  !fe(t,  •  )[(-,,  where  !l-||r  is  defined  in  (2.20),  we  note 
that  up  to  higher  order  terms  in  the  space  mesh  size, 

!!e(t ,  - )  ||£  -  !||e(t.  O'!!. 
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The  adaptive  procedure  used  in  the  experiments  is  that 
partly  described  in  section  1.  More  specifically,  it  is 
based  upon  controlling  :!  e(t,-)  for  all  t  [0,T-r,j..  ] 
according  to  the  following  principle.  For  any  given  tolerance 
TQL  >  0  and  constants  C1 ,  C„  satisfying 

(5.1)  0  <  C,  2  1/2  and  2  d  <  00 , 
the  space  meshes  should  be  modified  so  that 

(5.2)  !|e(t,*)  !|  ~  TOL;  It  6  [0,?_p.£T  ],  and 

(6.3)  Cx  •  TOL  5  iil  e  ( t ,  •  )  '  i  5  C^-TOL;  Iti  ^,-FINAL^- 

That  is,  it  is  desirable  that  l'|  e(  •  ,  •  )  ||  is  approximately 
equal  to  TOL,  but  acceptable  that  (6.3)  holds.  We  now 
make  a  number  of  remarks  concerning  underlying  assumptions 
used,  the  goals  (6.2)  and  (6.3),  and  the  computational  pro¬ 
cedures  employed  in  our  experiments  to  achieve  these  goals. 

(R.l)  The  adaptive  scheme  is  a  "once-through"  procedure. 

If  the  FEMOL  is  implemented  more  than  once  for  a  given 
problem,  then  no  information  obtained  from  a  previous  imple¬ 
mentation  is  used  adaptively  in  a  loop  to  solve  the  problem. 

(R.2)  The  primary  goal  is  to  control  '!  e ( -  ,  - )  j||  according 
to  (6.2)  and  (6.3),  and  not  to  minimize  it  for  some  fixed 
number  of  mesh  points,  as  is  often  proposed  in  adaptive 


schemes . 
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(R.3)  Decisions  concerning  mesh  modifications  are  based  upon 

achieving  (6.2)  and  (6.3)  with  |j;  e  (  -  ,  -  )  |  replaced  by  E  (  •  ) 

and  are  permitted  at  each  of  a  collection  {T  }  „  of 

m  m=0,d+l 

T 

times,  where  the  partition  1  =  { 0  =  TQ <  T^  < . . .<  TM+1  = 

TfiMAL ^  as  §^ven  a  Priori.  This  choice  of  possible  tran¬ 
sition  times  poses  no  actual  constraint  on  the  procedure,  and 
in  practice  these  points  of  time  could  include  output  target 
points  at  which  an  ODE  solver  must  be  restarted. 

(R.4)  The  choice  of  absolute,  and  not  relative  tolerance  in 
(6.2)  and  (6.3)  is  made  arbitrarily.  In  most  applications  the 
relative  tolerance,  or  a  combination  of  the  two  has  more  meaning. 

(R.5)  In  many  applications  it  may  be  acceptable  to  abandon 
(6.2)  and  place  a  low  priority  on  achieving  the  left  hand 
side  of  (6.3).  This  is  particularly  important  when  consider¬ 
ing  the  computational  work  involved  in  mesh  modification. 

(R.6)  TOL  is  assumed  to  be  time- independent .  In  applications 
involving  the  evolution  to  a  steady  state,  the  transient 
response  of  a  system  to  changing  loads,  or  the  propagation  of 
a  wave  through  different  media,  this  principle  can  be  modified. 

(R.7)  As  in  any  adaptive  procedure,  the  achievement  of  the 
goals  is  subject  to  certain  constraints.  In  our  experiments, 
the  number  of  constraints  was  kept  to  a  minimum  in  order  to 
more  easily  assess  the  potential  of  the  procedure.  The  primary 
constraints  were 
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(6.4)  M  <  N  ;  where  N  was  given  a  oriori, 

m  max  max  &  -  ’ 

(6.5)  h  <  H;  where  h  denotes  the  maximum  length  of  any 
finite  element  used  and  H  was  given  a  priori, 


and  in  modifying  a  partition  A  ,  only  combinations  of  mesh 
refinement  and  derefinement  were  permitted,  where 


(6.6) 


no  two  consecutive  nodal  points 
allowed  to  be  removed,  and 


m 

x .  _  , 

1--L 


x1?  were  both 
3 


(6.7) 


no  subinterval  ( x1?  .jX1?)  itself  was  permitted 

1-1  1  K 

refined  into  more  than  2U-  subintervals,  where 
positive  integer  q  was  given  a  priori. 


to  be 
the 


No  consideration,  however,  was  given  either  to  the  amount 
of  computational  work  needed  to  achieve  (6.2)  and  (6.3)  or  to 
the  question  of  how  best  to  balance  the  time  and  space  dis¬ 
cretization  errors.  Because  of  the  complex  interaction  of  the 
error  components  in  practical  computations  and  the  heuristic 
nature  of  many  work-related  constraints,  our  purpose  here  is 
to  isolate  the  space  discretization  error  and  examine  more 
clearly  the  potential  of  the  adaptive  procedure.  These  other 
important  issues  will  be  addressed  in  the  next  stage  of  the 
analysis.  The  ODE  per  step  tolerances  needed  to  essentially 
solve  the  resulting  systems  exactly  were  slightly  smaller 
than  those  used  in  the  experiments  in  [3]  with  time-independent 
meshes.  Because  of  this  stringent  requirement  and  the  orienta- 
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tion  of  our  experiments,  the  amount  of  CPU  time  used  in 
the  examples  is  irrelevant  and  will  not  be  quoted. 

Let  the  parameters  TOL,  C  ,  H,  Nmax>  and  q  be 

given.  For  integers  M  >  1  and  2  2,  suppose  that  the 

partitions  A x  =  { 0  =  TQ  <  T1  < .  .  .<  TM+1  =  TFIyAL'  and 
a  =  {0  =  xj?  <  x!?  < .  .  .<  x°  =1}  are  given,  and  that  the 

0  Ox  Nq 

variables  m  and  modflag  are  initially  set  equal  to  zero. 
The  algorithm  for  the  adaptive  construction  of  meshes  and 
the  FEMOL  solution  of  eqs .  (1.1)  can  then  be  summarized  as 


follows . 

ste; 

PA: 

(1) 

Determine  UCT^,*)  and  l 

(2) 

(3) 

Compute  {q  .}•  1  N  '  and 
J  J  ’  m 

If  the  exact  solution  u 

(4) 

wv 

Set  q° 

and  0  (T  )  . 

m 

1/2 

=  TOL/N*  ,  q  = 

ste 

£_1: 

If 

modflag  = 

1,  set  A*  =  A_ 
’  m  m 

A* 

m 

from  A  _ 
m 

by  checking 

(1)  If  j 

?  2-  'VV?1 

max{q 

m 

Xj-1 

and  marking  x? 

next 

pass . 

£  H,  and 
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(2)  If  ri  .  >  rj  ,  then  refine  (x^  ,  ,xrn)  into 
m3  1  -  3-1’  ] 

2^  uniform  subintervals,  where  integer 

p  6  {l,...,q}  approximately  solves  2^  =  n  ./q^. 

m3 


step  3: 


If  A"  i  Am  and  A”  defines  less  than  N  subintervals, 

mm  m  max  ’ 

then  set  Am  =  A”,  modflag  =  1,  and  redo  steps  1  and  2. 


stop;  otherwise  determine  U(t,*)  for 


t  €  (T  ,  T  ]  by  solving  an  N  -  1  dimensional  ODE 
m  m+1  °  m 

system. 


Set  A_. -i  =  A  ,  N  =  N  ,  modflag 
m+1  m  m+1  m  & 

steps  1-4  with  m  replaced  by  m+1. 


0;  and  redo 


In  step  1  of  each  loop  of  the  algorithm,  the  vector 

{U(T  jX1!1)}.  w  ..  of  nodal  values  of  the  FEMOL  solution  U 
m  3 

at  t  =  T  is  determined  as  the  solution  of  the  matrix  equa- 
tion  associated  with  the  standard  piecewise  linear  finite 
element  approximation  for  an  elliptic  boundary  value  problem. 
The  way  in  which  the  entries  of  these  matrices  and  the  indi¬ 
cators  {+1  M  were  computed  was  different  in  each  of 

m  j  i  “  j.  *  i  m 
J  J  m 

the  examples  and  will  be  discussed  shortly.  According  to  the 
manner  in  which  Ut(T  >•)  is  prescribed  (cf.  (2.35),  (2.36)) 
after  a  partition  Am  is  obtained  through  mesh  modification, 
the  computational  procedure  used  in  step  1  for  m  >  1  to 


obtain  1^(1^,  •  )  ,  and  hence  the  right  hand  side  of  the  matrix 
equation  is  the  following. 

In  the  event  that  mesh  derefinement  has  occurred,  i.e.g. 

a  point  x™  ^  has  been  removed  from  the  mesh,  a  system  of 

two  linear  equations  must  be  solved  to  determine  the  values 

U.(T  ,xra  )  and  U.(T  ,xI?)  on  the  new  mesh,  thus  determining 
t  m  j  -  2  t  m  3 

the  function  lh(T  ,■)  via  L„-proj ection .  If  mesh  refine- 
t  m  2 

ment  has  occurred,  simple  linear  interpolation  is  used  to 

determine  U.CT  ,  •).  While  these  processes  are  local,  the 
t  m 

process  for  obtaining.  U(T  ,  •)  is  not,  and  a  full  system  of 
N  -  1  linear  equations  must  be  solved.  The  new  mass  and 
stiffness  matrices  are  related  to  the  old  ones,  but  to 
simplify  the  data  structures  involved,  no  advantage  is  taken 
of  this  here  and  all  matrix  entries  are  recomputed  after  each 
mesh  modification. 


In  step  2  of  every  loop  of  the  algorithm  the  adaptive 

decision  making  process  takes  place.  Each  indicator  r)mj 

is  compared  to  the  root  mean  square  of  the  desired  error 

TOL.  The  underlying  assumption  used  here  is  that  the  best 

way  to  achieve  (6.2)  and  (6.3)  for  all  t  €  [T  ,T  , , )  is  by 
J  m  m+1  J 

constructing  A  so  that  the  indicators  {n  . } .  are 

m  ‘m3  3=l>Nm 

equilibriated .  This  assumption  is  based  upon  the  hypotheses 
which  led  to  Thms .  4.1  and  4.2,  and  is  made  in  view  of  results 
shown  in  [2],  where  optimal  meshes  for  elliptic  boundary  value 
problems  were  characterized  by  asymptotic  equality  of  the  error 
indicators . 


I 


The  mesh  modifications  implementing  this  principle  are 
restricted  to  being  combinations  of  refinement  and  derefine¬ 
ment  for  the  sake  of  computational  simplicity.  Since  each 

r)  .  is  proportional  to  |  x™  -  x1?  .  |  ,  it  is  assumed  that  by 
mi  i  1-1  y 

removing  the  nodal  point  x?  ^  in  constructing  A m+^  from 


A  we  will  obtain  q 
m  m+1  n 


0,  2  .2  v 1/ 2 

2  C  n  +  r|  .  ) 

m  j-1  mi 


where 


,  m+1  m+1*  ,,  m  m,  .  - . 

( Xn  i,xn  '  ~  ^xj  2,Xj  '  Similarly,  by  refining  any  interval 

,  m  m«  .  .  . , 

(x.  ,  ,x.)  into  the  union 
1-1  1 


vlt2P 


n-1  n 


of  2^  uniform  subintervals,  it  is  assumed  that 


vlt2P 


q  .  /  2^ . 
m] 


These  assumptions  are  based  upon  asymptotic  analysis. 

In  the  actual  implementation  of  step  2  of  the  algorithm, 

if  the  number  ^m+]_  subintervals  in  a  newly  constructed 

partition  A  ,,  differs  much  from  the  number  N  in  A  , 
m+ 1  mm 

then  each  of  the  previously  computed  indicators  {q  . ^ 
is  also  compared  to  the  values  of  q-  and  q  obtained  by 
replacing  N  by  N  , ,  .  This  is  done  to  avoid  a  subsequent 
mesh  modification  and  suggests  that  the  history  of  the  indica¬ 
tors  can  be  used  to  make  the  algorithm  more  efficient. 
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The  conducted  experiments  consisted  of  two  examples. 

In  all  runs  the  parameter  values  CL  =  1/2,  C_  =  2 ,  N  = 

1  2  3  max 

T 

200,  and  q  =  3  were  fixed.  The  partitions  and  A 

were  taken  to  be  uniform  and  defined  by 

(6.8)  x?  =  j/N  ;  j  =  0,  N  and 

3  0  0 

(6.9)  T  =  m  *  DT  =  m  *  (T.,,....  /M+l)  ;  m  =  0,  M  +  1. 

m  FINAL 

The  values  of  Ng ,  Tp^^,  M,  and  hence  also  DT  varied 
in  the  experiments,  as  did  the  parameters  TOL  and  H. 

Example  1. 

The  coefficients  a  and  b  and  the  exact  solution  u 
of  eqs.  (1.1)  are  as  in  Example  1  of  [3,  section  7],  namely 
a(x)  =  cosh(4x-2),  b(x)  =  sinh(2x),  and 

(6.10)  u(t,x)  =  u1(t,x)  +  [1/2  +  1/2  tanh( 10t-ll) ]u^ ( t , x) ,  where 

(6.11)  u1(t,x)  =  1  -  exp  {-x"V(10t+  .1)}  and 

2 

(6.12)  u  (t,x)  =  2  sin(TTx)  +  2  cos (  2n t ) sin(  2nx)  . 

For  0  s  t  2  .9,  u1  dominates  the  behavior  of  u  and 

u  decays  in  t.  After  undergoing  a  sharp  transition  for 

.  .  .  2 
.9  <  t  <  1.1,  u  oscillates  as  its  dominant  component  u 

(cf .  figure  1) . 

In  this  example  standard  two-point  Gaussian  quadrature 
was  used  in  the  computation  of  the  mass  and  stiffness  matrix 
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entries,  load  vectors,  error  indicators,  and  absolute  error 
Etru-  The  parameter  values  NQ  =  10,  TFINAL  =  3, 

M  =  49(=»  DT  =  .06),  and  H  =  .2  were  fixed,  and  the  adaptive 
FEMOL  procedure  was  implemented  twice,  with  the  absolute  error 
tolerance  TOL  taken  to  be  .530  and  .265.  These  are  the 


mean  values  about  which  En 


oscillated  for  1  <  t  <  3 


Example  1  of  [3,  section  7],  using  20  and  40  uniform,  t- 
independent  finite  elements  (cf.  [3,  figure  4],  where  the 
corresponding  relative  error  ERRF  is  shown  oscillating  about 
.05  and  .025). 


The  numbers  of  elements  used  as  functions  of  time  for  the 


different  values  of  TOL  are  shown  in  figure  2.  Because  no 
constraint  prohibiting  minor  mesh  modifications  was  imposed, 
mesh  changes  occurred  at  most  of  the  50  allowed  transition 
times . 


Figure  3  illustrates  the  effectivity  of  the  adaptive 
procedure  in  this  example.  For  each  TOL 


(6.13)  1/2  ■  TOL  2  ETRU(t)  <  2  •  TOL;  I  <  t  <  3. 

For  0  <  t  <  1,  the  imposition  of  the  maximum  mesh  size  con¬ 
straint  limited  the  flexibility  of  the  procedure,  and  there¬ 
fore  ETR1J(t)  «  TOL. 

The  effectivity  of  the  estimator  E(-)  for  meshes 
changing  according  to  the  adaptive  procedure  is  illustrated 
in  figure  4.  For  each  TOL 


p 
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(6 

.14) 

.95 

<  0(t)  £  1.12; 

Vt  € 

[.9,1 

.13, 

(  6 

.15) 

.99 

<  0(t)  <  1.06; 

Vt  € 

[1.1, 

33, 

it  appears  as  though  0 Ct)  -*■  1  as  TOL  -*■  0  for  most 
t  €  [1.1,3],  The  fact  that  0(t)  is  not  as  near  1  for 
t  €  (0,.9)  is  probably  due  to  the  flatness  of  the  exact 
solution  u(t,*),  which  essentially  violates  an  assumption 
in  the  theory  of  section  4. 

Example  2 . 

The  exact  solution  u  of  eqs .  (1.1)  is  as  in  Example  2 
of  [3,  section  7],  namely 

(6.16)  u(t,x)  =  1/2  +  1/2  tanh[2p(x  -  lOt) J;  ?  -  20. 

The  solution  is  a  wave  with  approximate  front  width 
(3  ^  which  moves  in  the  positive  x-direction  at  speed  10, 
and  is  as  pictured  in  figures  5(a),  (b),  (c). 

The  coefficients  a  and  b  of  eqs.  (1.1)  were  taken 
here  to  be  the  constant  functions  a  =  1  and  b  =  0.  This 
choice  was  made  so  that  all  integrations  needed  in  the  compu¬ 
tation  of  the  matrix  entries,  load  vectors,  error  indicators, 
and  could  be  performed  exactly.  Strong  pollution 

effects  were  observed  in  the  adaptive  procedure  for  this 
example  when  two-point  Gaussian  quadrature  was  employed. 

This  effect  was  particularly  evident  in  computing  the  L^- 
projections  of  U  needed  to  determine  the  initial  data  U 
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at  the  transition  times,  and  indicates  the  care  with  which 
numerical  integration  methods  should  be  chosen  when  applying 
the  adaptive  FEMOL  procedure. 


The  parameter  was  se1:  eaua-*-  to  -04>  which  is 

the  amount  of  time  needed  for  the  center  of  the  wave  to 


travel  from  x  =  0  to  x  =  .4.  Two  implementations  of  the 
adaptive  scheme,  which  we  shall  refer  to  as  Examples  2(a)  and 
2(b)  were  carried  out.  In  Example  2(a),  the  remaining  para¬ 
meters  were  taken  to  be  NQ  =  10,  M  =  19(=  DT  =  .002), 

H  =  .1,  and  TOL  =  .8.  The  values  taken  in  Example  2(b) 
were  NQ  =20,  M  =  39 (  DT  =  .001),  H  =  .05,  and  TOL  =  .1. 
The  desired  absolute  errors  of  .8  and  .1  correspond 
approximately  to  16%  and  2%  errors. 

In  figure  6(a)  the  nodal  points  at  every  fifth  transition 
time  in  the  first  half  of  the  interval  (0,1)  chosen  by  the 
adaptive  scheme  in  Example  2(a)  are  pictured.  The  total 
number  of  nodal  points  at  any  time  varied  between  16  and 
18.  In  figure  6(b)  the  nodal  points  at  every  tenth  transition 
time  for  Example  2(b)  are  shown.  Between  64  and  70  nodal 
points  were  used  to  partition  (0,1)  at  each  time  in  this 
example . 

As  in  Example  1,  in  both  Examples  2(a)  and  2(b)  a  mesh 
change  occurred  at  nearly  every  allowed  transition  time. 
Figures  6(a)  and  6(b)  clearly  illustrate  the  ability  of  the 
adaptive  scheme  to  concentrate  nodal  points  along  the 
characteristic  x  =  lot  of  the  wave,  where  all  of  the 
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activity  takes  place. 

The  effectivity  of  the  adaptive  procedure  is  illustrated 
in  figure  7.  Here,  as  also  in  figure  3,  the  output  time 
points  were  chosen  to  coincide  with  the  allowed  transition 
time  points.  After  the  second  mesh  transition  in  Example 


Etru  oscillates  about  the  value  .55  and  we  see  that 


(6.17)  1/2  •  TOL  £  ETRU(t)  £  TOL. 


Similarly,  ETPI1  oscillates  about  .08  in  Example  2(b) 


after  the  second  mesh  transition  and 


(6.18)  3/4  TQL  5  ETRU(t)  £  5/4  T0L- 

The  effectivity  of  the  estimator  E(-)  in  both  Examples 
2(a)  and  2(b)  is  illustrated  in  figure  8.  In  Example  2(a) 
at  all  but  one  of  the  21  time  points,  0  satisfies 

(6.19)  1  £  0(t)  £  1.1, 

and  at  all  of  the  41  time  points  in  Example  2(b) 

(6.20)  .95  £  0(t)  £  1.05. 

The  steady  decline  of  0  as  a  function  of  t  in  Example 

2(b)  is  apparently  due  to  errors  introduced  by  the  discontin¬ 

uous  mesh  transitions.  In  repeating  Example  2(b)  with 
H  =  .1  instead  of  .05,  we  observed  a  similar  decline  in  the 
effectivity  ratio  0,  but  this  time  to  .85.  In  another 
experiment  with  H  =  .1  and  DT  =  .002  instead  of  .001, 
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this  decline  was  slightly  alleviated. 

These  experiments  suggest  modifying  the  estimator  E(0, 
as  was  proposed  in  section  5  (cf.  (5.45)),  when  no  prior 
knowledge  of  the  relationship  between  TOL  and  H,  or  TQL 
and  DT  is  available.  We  recall  that  the  theory  in  section 
4  concerning  the  effectivity  of  E(-)  was  based  partly  upon 
assumptions  related  to  the  ratios  of  the  sizes  of  elements 
in  each  mesh  and  the  distances  between  mesh  transitions. 

As  was  observed  in  [3]  for  t-independent  meshes,  it 
has  been  seen  in  the  selected  examples  presented  here  that 
E(*)  is  an  effective  estimator  of  E„OIt.  While  effectivity 
ratios  anywhere  in  the  range  [.1,10]  would  be  quite 
acceptable  in  most  applications,  0  was  very  near  one  for 
all  of  the  adaptively  constructed  meshes.  In  view  of  the 
apparently  stringent  conditions  concerning  the  asymptotic 
range  for  mesh  sizes,  which  were  given  for  general  changing 
meshes  in  the  theory  of  section  4,  this  is  a  particularly 
encouraging  result. 

The  mesh  modification  scheme  presented,  in  which  decisions 
based  upon  E(*)  are  made  adaptively  by  the  computer,  was  like¬ 
wise  quite  effective  in  our  experiments.  In  each  case  the  goals 
of  the  procedure  were  realized,  thus  illustrating  the  robustness 
of  both  E(*)  and  the  adaptive  process.  These  results  from  the 
selected  examples  suggest  the  method's  applicability  to  more 
general  problems  occurring  in  practice  (cf.  [4]). 
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Figure  1 


X 

Exact  Solution  u(T,x)  of  E 
X  for  various  T . 
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Figure  4.  Effectivity  Ratio  Q(T)  in  Example  1  vs .  T 
TOL  =  Absolute  Error  Tolerance  . 
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Figure  5(a),(b),(c). 


Exact  Solution  u(T,X)  of  Example  2 
vs.  X  for  various  T. 


Figure  6(a).  Nodal  Points  in 
Every  5th  Grid  in  Example 
2(a)  . 
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Figure  6(b).  Nodal  Foints  in 
Every  10th  Grid  in 
Example  2(b) . 
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